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A  FTNITK  ELEMENT  ANALYSIS  OF  SHOCK  AND  FINITE-EMPLITUDE  WAVES 
IN  ONE -DIMENSIONAL  HYPERELASTIC  BODIES  AT  FINITE  STRAIN 

R.  B.  Fost  and  J.  T.  Oden 

Department  of  Engineering  Mechanics,  The  University  of  Alabama  in  Huntsville 

SUMMARY.  The  general  theory  of  shock  and  acceleration  wc 'es  in  isotropic, 
incompressible,  hyperelastic  solids  is  used  in  conjunction  with  the 
concept  of  finite  elements  to  construct  discrete  models  of  highly  non¬ 
linear  wave  phenomena  in  elastic  rods.  A  numerical  integration  scheme 
which  comb i nes  features  of  finite  elements  and  the  Lax-Wendroff  method 
is  introduced.  Numerical  calculations  of  the  critical  time  for  shock 
formulation  are  given.  Numerical  results  obtained  from  representative 
cases  are  discussed. 

1.  INTRODUCTION 


Until  very  recent  times,  the  quantitative  study  of  the  dynamic 
response  of  highly  elastic  solids  at  finite  strain  has  stood  outside  the 
reach  of  existing  analytical  and  computational  methods.  The  complete 
dynamical  theory  is,  itself,  still  being  pieced  together  some  20  years 
after  the  modern  era  of  finite  elasticity  began  (e.g.  [1,2]),  and  while  a 
few  solutions  to  finite-amplitude  vibration  problems  have  been  contributed 
(e.g.  [3, A, 5]  or,  for  additional  references  [oj),  and  while  some  advances 
have  been  made  in  the  theory  oi  finite  elastic  waves  (e.g.,  [7,8,9]),  actual 
calculations  invariably  involve  rather  ideal  geometries,  boundary-  and 
initial  conditions,  and/or  material  properties.  The  highly  nonlinear 
character  of  the  m'H.ient'.m  equatiers  ft  *  the  most  simple  hypere.lastic 
material  does  not  acco-nt  for  all  of  t  e  co.  putatione.  problems--by  definition. 


A 
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the  hyperelastic  solid  possesses  no  dissipative  mechanism  to  provide 
snoothing  or  damping  of  higher  frequencies.  Consequently,  the  computa¬ 
tionally  convenient  features  of  damping  enountered  in  nonlinear  visco¬ 
elasticity  and  thermoviscoelasticity  calculations  [10,11,12]  are  not 
present.  To  complicate  matters,  it  is  now  generally  recognized  that  shock 
waves  can  be  easily  produced  in  such  materials,  even  when  smooth  initial 
conditions  are  prescribed.  The  recent  experimental  work  of  Kolsky  [13] 
gives  evidence  to  the  possibility  of  even  tensile  shock  waves  developing 
in  certain  stretched  natural  rubbers,  a  phenomena  already  anticipated  in 
the  theoretical  work  of  Bland  [8]  end  Chu  [9].  Tn  such  cases,  practically 
all  of  the  popular  numerical  integration  schemes  now  used  in  structural 
dynamics  are  ineffective. 

In  the  present  paper,  we  consider  an  important  subclass  of  the  prob¬ 
lems  alluded  to  in  our  opening  comments,  and  we  demonstrate  that  the  finite- 
element  method,  when  used  in  conjunction  with  established  methods  of  com¬ 
putational  hydrodynamics  and  implemented  with  modern  computing  machinery, 
can  he  extremely  effective  for  this  class  of  problems.  More  specifically, 
tin-  present  paper  is  concerned  with  the  application  of  the  finite-element 
method  to  the  calculation  of  shock  waves  and  finite-amplitude  acceleration 
waves  in  highlv  elastic  rubber-1  Ike  materials.  So  as  not  to  obscure  con¬ 
ceptual  and  physical  details,  this  investigation  is  confined  to  the  study 
of  only  longitudinal  motions  of  finite  homogeneous  rods  of  Isotropic,  In¬ 
compressible,  hyperelast i c  materials.  We  hope  to  treat  the  more  difficult 
problems  of  finite-amplitude  waves  in  two-  and  three-dimensional  hyperelos- 
tic  bodies  in  later  work.  Herein,  we  develop  a  fully  discrete  representa¬ 
tion  of  the  behavior  which  employs  a  finite-element  approximation  of  the 
spatial  variation  of  the  displacement  field  or  the  longitudinal  extension 
ratio  and  a  finite-difference  representation  of  the  temporal  behavior.  For 


shock  cal culat ions ,  we  introduce  an  apparently  new  explicit  integration 
scheme  which  is  shown  to  he  very  effective  in  handling  the  formation, 
propagation,  and  reflection  of  shocks  in  disturbed  media.  The  scheme 
involves  t he  use  of  finite-elements  as  a  basis  for  formulating  I.ax-Wendrof  f- 
type  time  integration  algorithms.  Moreover,  we  also  consider  the  problem 
of  computing  numerically  various  critical  times  required  for  the  formation 
of  shocks  from  smooth  and  Lipschitz-cont inuous  initial  data.  In  the  final 
section  of  the  paper,  we  cite  numerical  results  obtained  from  a  number  of 
representative  test  cases. 

2.  PHVSTCS  OF  WAV1-S  IN  NONLINEAR  K  LA  STIC  MATKRIALS 

In  this  section,  we  shall  review  certain  features  of  the  dynamical  and 
thermodynamical  theory  of  finite  elasticity  that  are  essential  to  our 
study.  More  complete  details  can  be  found  in  the  books  of  Green  and  Zerna 
[14],  Green  and  Adkins  (.15],  and  Bland  [8],  and  certain  special  features 
are  discussed  in  the  papers  of  Nowinski  r 16 3 ,  Ames  [17],  and  Reddy  and 
Achenbach  [18]. 

2.1  Motion  of  a  Thin  Rod.  We  begin  by  considering  longitudinal  motions 
of  a  thin  rod  of  isotropic,  incompressible  material.  While  at  rest  in  a 
natural  reference  configuration  rto ,  the  rod  has  a  uniform  symmetric  cross 
section  of  area  Aq  ,  a  length  T,0 ,  and  a  mass  density  Po.  To  describe  the 
motion  of  the  rod  relative  to  its  reference  configuration,  we  establish  a 
fixed  spatial  frame  of  reference  x,  with  origin  0  at  one  end  of  the  bar, 
which  is  assumed  to  be  restrained  from  motion.  We  denote  by  X  the  labels 
of  material  particles  (material  coordinates)  of  the  bar,  and  we  select 
these  labels  so  as  to  numerically  coincide  with  the  spatial  coordinates  x 
when  the  bar  occupies  its  reference  configuration.  The  function  x  =  x(X,t) 
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then  gives  the  spatial  position  of  the  particle  X  at  time  t  and  defines 
the  longitudinal  motion  of  the  bar,  while  u(X,t)  =  x(X,t)  -  X  defines  the 
displacement  of  particle  X  at  time  t. 

The  material  gradient  dx/dX  serves  as  a  convenient  measure  of  the 
deformation;  physically,  it  represents  the  longitudinal  extension  ratio  X 
(also  called  the  stretch)  which  is  expressible  in  terms  of  the  displacement 
gradient  according  to 

x  -1  = 1  +  !  ■ 1  +  “■  <2-'> 

where  here,  and  henceforth,  we  use  the  subscript  notation  ux  =  du/dX  to 
denote  partial  differentiation  with  respect  to  X.  The  exlensional  strain 
Y  is  then  simply  (X2  -  l)/2. 

Now  any  disturbance  supplied  to  the  rod  will  travel  with  finite  speed 
from  one  particle  to  another  in  the  form  of  a  surface  (curve)  of  discon¬ 
tinuity  s  (a  wave)  in  the  X-t  plane.  If  we  denote  by  Y(t)  the  particle 
X  reached  by  a  wave  front  at  time  t,  then  the  set  of  points  (Y(t),t),  r  being 
a  real  parameter,  defines  a  cruve  in  the  X—  t  plane  across  which  jump 
discontinuities  in  various  partial  derivatives  of  u(X,t)  can  occur.  Indeed, 
acceleration  waves  involve  jumps  in  the  acceleration  ult  (and  the  stress 
gradient  do/dX)  and  shock  waves  (shocks)  involve  the  propagation  of  sur¬ 
faces  across  which  jump  discontinuities  in  the  velocity  u(X,t)  =  ut  =  du/dt 
and  the  stress  are  experienced.  The  intrinsic  wave  speed  relative  to  the 
material,  denoted  c,  is  then  simply  dY(t)/dt.  However,  the  spatial  position 
y(t)  of  the  wave  is  clearly 

y(t)  =  x(Y(t) , t)  (2.2) 

and  the  absolute  wave  speed,  as  seen  by  a  stationary  observer,  is 

v=  _  dx(Y(t).t) 

dt  dt 


(2.3) 


5 


2.2  The  Balance  l,<iws.  We  assume,  of  course,  that  the  response  of  bar 
satisfies  the  basic  conservation  laws  of  physics.  For  the  problems  at 
hand,  these  assume  the  following  global  forms: 


etc.  Mass  is  conserved  in  the  rod  and  we  have  ignored  body  forces  and 
internal  heat  sources  for  simplicity. 

At  particles  X  that  do  not  fall  on  a  surface  of  discontinuity,  (2.4)- 
(2.6)  lead  to  the  local  forms  of  the  balance  laws 


() 

P0 U  -  o,  =  0 

c  -  oX  -  qx  =  0  (2.8) 

-  q*  +  —  o 

whereas  at  the  surface  of  discontinuity,  we  obtain  the  jump  conditions 

p0c£u]!y  4  £oHy  =  o 

\  PBclV]y  +  c  £<-])  Y  +  (L CJU u  Y  +  I£qH  Y  =  o  (2.9) 

ctSHY  -  S!g]!y  •'  0 

It  is  often  convenient  to  introduce  the  Helmholtz  free  energy  cp(X,t) 
and  the  internal  dissipation  6(X,t)  defined  for  X  t  Y(t)  by 

cp  =  e  -  £8  and  6  =  8£  -  qx  (2.10) 

Then  the  last  member  of  (2.8)  can  be  rewritten  in  the  alternate  form 

i0.  ...  i0. 

In  addition,  we  can  also  impose  Maxwell's  theorem,  which  asserts 
that  for  any  function  f(X,t)  jointly  continuous  in  X  and  t,  but  whose  first 
partial  derivatives  fx  and  f  suffer  jump  discontinuities  at  S,  the  jump 
across  S  in  the  (two-dimensional)  gradient  of  f  must  be  parallel  to  a 
vector  (-1,  c)  normal  to  S.  Applying  this  idea  to  the  motion  x(X,t) 
yields  Hadamard's  compatibility  condition: 

luly  +  c{[\3y  =  0  (2.12) 


2.3  Waves  in  Hyperelastic  Materials.  We  now  aim  our  analysis  toward 


waves  in  hyperelastic  materials;  that  is,  we  wish  to  consider  materials 
for  which  there  exists  a  potential  W  which  is  a  function  of  the  current 
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value  of  X,  and  for  which 

•  •  dw 

W  =  oX  and  a  =  |?  (2.13) 

OA. 

The  question  arises , however,  as  to  whether  or  not  a  theory  of  hyperelasticity 
is  reconcilable  within  the  thermodynamic  framework  established  thusfar. 

This  is  a  classic  question,  and  standard  arguments  can  be  found  in  a  num¬ 
ber  of  places  (e.g.  [6,12])  to  the  effect  i.  >  .yperelasticity  is  Indeed 
possible  in  a  num^r  of  physically  meaningful  situations.  The  fact  that 
these  standard  arguments  are  not  valid  at  surfaces  of  discontinuity  is 
fundamental  to  the  physics  of  shock  waves. 

We  mention  two  cases.  First,  consider  a  class  of  perfect  materials 
(cf. [6],  p.  296),  the  constitution  of  which  is  defined  by  equations  for  e, 
cr,  0,  and  q  depicted  as  functions  of  the  current  values  of  X  and  5,  with  q 
also  dependent  on  9x(X,t).  For  reversible  processes  performed  on  such 
materials,  the  dissipation  6  =  9§  -  qx  *  0,  and  (2.8)a  gives 

(°  -  S) u  (e  -  It)  5 + 1  «9. 2  0  <2-14> 

so  long  as  X  e  (Lc  -  Y(t)).  Secondly,  we  consider  a  class  of  simple 
materials  (cf.  [12],  p.  202)  whose  constitution  is  defined  by  equations  for 
cp,  o,  §,  and  q  in  terms  of  current  values  of  X  and  9,  with  q  dependent 
upon  9x(X,t)  also.  For  reversible  process  (6  =  0),  using  (2.10)  in  (2.8) 
gives  us  the  inequality 

for  all  X  €  (L0  -  Y(t)).  If  these  two  inequalities,  (2.14)  and  (2.15),  are 
to  be  maintained  for  arbitrary  rates,  it  necessarily  follows  (cf  (12],  p. 
214)  that  in  the  absence  of  a  shock,  a  theory  of  hyperelasticity  is 
appropriate  for  reversible  isentropic  processes  (6=0,  §  =  0)  performed  on 
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the  above  class  of  perfect  materials  and  for  reversible  isothermal  processes 
(b  -  0,  3  =  0)  performed  on  the  above  class  of  simple  materials.  Tn  the 
former  case,  the  strain  energy  is  associated  with  the  internal  energy,  in 
the  latter  ease  it  is  associ  ted  with  tiie  free  energy.  However,  since  a 
s».  ek  is  characterized  by  discontinuities  in  the  displacement  gradient, 

‘tie  necessary  derivatives  of  e  in  (2.14),  or  of  <f>  in  (2.15),  do  not  exist 
for  X  =  Y(t).  Therefore,  we  must  have  energy  dissipation  at  X  =  Y(t),  i.e., 
6  =  95  -  qx  >  0,  and  we  lose  the  notion  of  reversibility. 

Due  to  the  considerable  difficulties  involved  in  solving  the  nonlinear 
thermo-mechanical  equations  governing  irreversible  thermodynamic  processes, 
the  only  exact  solutions  available  (cf.  [19,20])  are  for  shocks  with  uniform 
conditions  on  both  sides  of  the  oiscontinuity .  Hence,  for  additional 
solutions,  we  need  to  simplify  the  governing  equations  so  that  they  become 
tractable.  One  possibility,  suggested  by  the  exact  discontinuous  solutions 
themselves,  is  the  well  known  fact  that  "Weak  shocks"  are  nearly  isentropic 
(e.g.,  see  [8],  [2l],  or  [22]  for  detailed  discussions).  That  is,  taking 
the  proportional  change  in  magnitude  across  the  shock  of  some  state  para¬ 
meter,  say  ux,  as  a  measure  of  shock  "strength",  the  change  in  entropy 
across  the  shock  is  only  of  third  order  in  the  shock  strength  for  small 
changes  of  ux.  Therefore,  for  weak  shocks,  we  can  neglect  the  entropy 
change  and  consider  5  as  constant  for  all  X  and  t,  i.e.,  the  deformation 
takes  place  isentropical ly.  It  is  of  special  interest  to  consider  this 
"isentropic  approximation"  when  the  initial,  or  reference  configuration 
is  the  natural  stress-free  state  where  X(X,0)  =  1,  9(X,0)  =  T0  =  constant. 
Then  the  isentropic  approximation  becomes  5=0  everywhere  for  all  time. 
Moreover,  for  this  case,  the  isentropic  approximations  render  the  mechsuLcal 
equations  independent  of  the  thermal  equations,  and  the  mechanical  jump 
conditions  (the  first  of  (  2,9)  and  (  2.12))  alone  are  sufficfent  to 
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determine  the  shock  process.  (Naturally  the  energy  jump  condition  remains 
valid,  but  here  it  would  only  be  used  to  check  the  energy  balance  after 
solving  the  problem.)  Also,  we  can  readily  define  the  strain  energy 
function  W  of  (2.13)  in  terms  of  the  internal  energy  W(\)  *  e(^»5)l^mQ*  8° 
that  we  have  the  constitutive  relation  for  the  stress  c  =  de/iX.,  in  agree¬ 
ment  with  (2.13). 

We  'eroark  that  the  local  balance  laws  (2.8)  suggest  another  simplifica¬ 
tion:  the  adiabatic  approximation,  wherein  ve  assume  the  heat  conduction 
small  enough  to  take  q  «  0.  Outside  the  shock  region,  the  adiabatic  pro¬ 
cess  is  reversible  (6  *  0);  and  then  from  (2.10)  we  get  the  reversible 
adiabatic  process  to  be  an  isentroplc  process.  Hence,  in  regions  of  the 
rod  where  X  ^  Y(t),  a  theory  of  hyperelasticity,  in  the  sense  of  (2. 14),  is 
possible.  At  the  shock,  the  adiabatic  process  is  not  reversible:  entropy 
is  produced.  Then,  by  integrating  the  inequality  in  (2.8)  with  q  *  0,  we 
get  §  -  5(X).  Therefore  the  entropy  at  each  material  point  has  a  constant 
value  unless  a  shock  passes  over  the  point,  at  which  time  the  value  of  the 
entropy  is  changed  to  a  new  constant.  Consequently,  until  the  time  at 
which  a  shock  forms,  the  adiabatic  process  is  isentroplc  for  every  X  e  (0,Lo); 
after  this  time,  the  adiabatic  process  is,  in  general,  piecewise  isentroplc, 
i.e.,  it  is  isentroplc  for  every  X  e  [(0,Y(t")]»  [Y(t+),Lo)}. 

Both  of  these  assumptions  can,  «»f  course,  be  simultaneously  incorporated 
if  we  consider  the  propagation  of  weak  adiabatic  shocks.  Tn  this  case, 
again  following  (2.14),  the  material  will  be  everywhere  hyperelnstic  at 
all  times. 

3,  EVOLUTION  AND  PROPAGATION  OF  DISCONTINUITIES 


In  this  section  we  will  look  at  the  physical  conditions  which  are 
generally  required  for  the  formation  and  propagation  of  discontinuities — 

both  shock  waves  and  acceleration  (simple) waves  in  rubber-like  materials. 
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We  also  briefly  consider  methods  for  determining  the  time  required  for 
discontinuities  to  evolve  during  the  solution  process. 


3.1  Propagation  of  Shock  and  Acceleration  Waves.  In  the  absence  of  shocks, 
we  previously  obtained  the  first  member  of  (2.8)  as  the  local  form  of  the 
law  of  conservation  of  linear  momentum.  In  view  of  the  constitutive 
relation  (2.13)  for  hyperelastic  materials  we  have  o  =  0( X) #  so  that  the 

local  momentum  equation  for  such  materials  can  be  written  in  the  form 

i 

u  -  ca(ux)utx  =  0  (3.1) 


where  the  squared  intrinsic  wave  speed,  c2(ux),  is  given  by 


2 /  \  1  do 

c2(u  )  =  — 

*  Po  QX 


(3.2) 


We  also  note  that,  since  X  =  ux,  we  can  recast  the  local  momentum  equation 
(3.1)  in  terms  of  X  according  to 


*X  =  lc*(X)Xjx  (3.3) 

where,  clearly,  the  forms  of  ca(X)  and  c2(ux)  will  be  different. 

As  noted  earlier,  for  hyperelastic  solids  the  stress  is  derivable 
from  a  potential  function  W  which  represents  the  strain  energy  per  unit 
undeformed  (reference) volume  Uo .  For  isotropic  incompressible  bodies,  W 
is  generally  defined  in  terms  of  the  first  two  principal  invariants,  I1  and 
Is ,  of  Green's  deformation  tensor,  the  third  principal  invariant  being 
unity.  In  the  present  case,  =  2X  +  X  ,  I2  *  2X  +  X"2,  and  elimination 
of  the  hydrostatic  pressure  with  the  condition  that  transverse  normal 
stresses  are  zero,  leads  to  the  general  constitutive  law 


o  =  2 (V.'j  X  +  WE)(1  -  X-3)  (3. A) 

where  W^  =  dw/dl^,  or  »  1,2.  Substituting  (3.4)  into  (3.2),  we  observe  that 
the  square  of  the  wave  speed  in  materials  defined  by  (3. A)  is  of  the  form 
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c3  =  J-  [(1  +  2X'3)W1  +  3X“*W?  +  2(1  -  X“*)8(Wll18  4  2W1P\4  W„)] 

*0 

(3.5) 

If  we  eliminate,  on  physical  grounds,  the  possibility  of  complex 
wave  speeds,  we,  in  turn,  impose  conditions  on  the  form  of  W  and  its  deriva¬ 
tives  W^,  W^j.  In  this  regard,  we  shall  assume  that  the  stiess  0  is  a 
continuous  monotonieally  increasing  function  of  the  stretch  X,  so  that  for 
all  X  €  (0,“),  we  have  0  <  p0c3(X)  =  da/dX  <  This  important  property 
allows  us  to  interpret  qualitatively  a  number  of  interesting  nonlinear  wave 
phenomena.  For  example,  suppose  that  a  time-dependent  surface  traction  is 
applied  at  the  free  end  of  the  rod.  During  each  infinitesimal  increment  in 
time,  the  corresponding  increment  in  load  produces  a  "wavelet,"  so  that, 
using  Nowinski's  terminology  [16],  the  net  effect  of  the  loading  is  to 
produce  an  "infinite  manifold"  of  wavelets  propagating  along  the  rod. 
Obviously,  each  successive  wavelet  propagates  at  a  speed  determined  by 
the  instantaneous  slope  of  the  0  vs.  X  curve  for  the  material.  Thus  if 
consecutive  wavelets  are  propagated  with  decreasing  speeds,  the  slope  of 
the  wave  front  will  gradually  decrease  (contrasting  markedly  with  the 
usual  sharp  discontinuity  at  the  wave  front  in  materials  with  linear  o  -  X 
curves),  and  the  response  will  be  propagated  as  a  simple  wave.  However, 
if  the  distance  between  successive  wavelets  decreases  during  propagation 
(they  are  generated  with  increasing  speeds),  the  wave  profile  steepens 
until  the  discontinuity  is  transformed  into  a  shock  wave. 

To  be  more  specific,  consider,  for  example,  the  following  special 
forms  of  the  strain  energy  function: 

(i)  The  neo-Hookean  form,  W  =  C(It  -  3)  (3.6) 

(ii)  The  Mooney  form,  W  =  C1(Il  -  3)  +  -  3)  (3.7) 

(iii)  The  Biderman  form,  W  =  B; (1^  -  3)  +  B2  (Ij  -  3)2 

+  Ba(Il  -  3)3+B4(I2  -  3) 


i 


i 

i 


i 

l 

! 


i 


( 

t 

I 

i 

! 


4 

J 


(3.8) 
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Here  C,  C1 ,  C2,...,B4  are  material  constants.  Examples  of  a  variety  of 
other  forms  of  W  proposed  for  real  materials  are  summarized  in  [12],  Note 
that  for  Mooney  materials 

c2  »  —  [Cx(l  +  2X"3)  +  3Cj, X"4]  (3.9) 

P0 

whereas  in  the  case  of  Riderman  materials 

c2  =  —  {[Bl  +  2R?(2X“1  +  X2  -  3)+  3Rj (2X-1  +  X2  -  3)s](l4  2X'3) 

Po 

+  3B4X-4  +  [4b2  +  12B3(2X-1  +  Xs  -  3)](X  -  X-2)2}  (3.10) 

The  function  c 7  for  neo-Hookean  materials  follows  from  (3.9)  by  setting 
Cj,  =  0.  Equations  (3.6)-(3.8)  with  (3.4)  and  (3.9)  and  (3. 10)  describe 
a  -  X  and  c  -  X  curves  of  the  type  shown  in  Figs.  1  and  2.  Clearly,  the 
type  of  wave  generated  by  an  initial  disturbance  depends  upon  both  the 
initial  state  (i.e.,  the  initial  value  of  X)  and  whether  X  is  subsequently 
increased  or  decreased.  A  discontinuity  is  propagated  as  a  simple  wave 
if,  and  only  if,  the  intrinsic  wave  speed  of  the  material  in  front  of  the 
discontinuity  is  greater  than  that  of  the  material  behind  the  discontinuity 
(cf.  [21],  p.  243).  Hence,  Fig.  2  suggests  that  for  both  the  Mooney  and 
neo-Hookean  materials,  only  compression  shock  waves  '-an  be  developed. 
However,  for  certain  Biderman-type  materials  (curve  (',  ‘n  Fig.  2),  it  is 
possible  to  produce  a  tensile  shock  wave  if  the  material  ia  preatretched 
sufficiently.  The  development  of  such  tensile  shocks  has,  in  fact,  been 
observed  experimentally  by  Kolaky  t 1 3 J .  It  is  also  interesting  to  consider 
the  case  in  which  an  applied  tension  load  is  suddenly  removed.  From  Fig. 

2  we  can  see  that  this  discontinuity  cannot  be  propagated  by  a  simple  wave 
since  the  wave  speed  is  smaller  before  the  load  is  released  than  after. 
Thus,  it  follows  that  the  instant  the  stress  at  the  end  of  the  rod  begins 
to  decrease  (after  having  increased)  is  also  the  time  at  which  the  first 


/ 


Fig.  2.  Wave  speed  versus  stretch  for  some  rubber-like  materials. 
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wavelet  emanates  from  the  end  with  a  propagation  speed  faster  than  the 
preceding  wavelets.  Therefore,  at  some  time  subsequent  to  the  moment  when 
the  unloading  starts,  the  propagated  influence  of  this  release,  or  un¬ 
loading,  is  expected  to  develop  into  a  shock  wave. 

3.2  Evolution  of  Discontinuities.  The  formation  of  discontinuities  in 
solutions  of  nonlinear  hyperbolic  equations  has  been  a  topic  of  interest 
for  many  years.  The  research  on  this  topic  falls  into  two  principal  cate¬ 
gories:  discontinuities  which  evolve  from  Lipschitz  continuous  initial 
data  and  those  which  evolve  from  smooth  initial  data.  In  the  former  case, 
there  is  a  clearly  defined  wave  front  and  a  characteristic  along  which  the 
initial  discontinuity  (even  when  starting  with  analytic  initial  data,  a 
Lipschitz  discontinuity  can  subsequently  develop)  propagates  until  it  tends 
to  a  Jump  discontinuity  at  some  critical  time,  say  tc„ .  This  jump  discon¬ 
tinuity  then  propagates  in  a  completely  different  manner  from  the  Lipschitz 
discontinuity.  As  this  is  discussed  in  detail  in  [8],  [21],  and  [23],  we 
will  mention  only  the  essential  features. 

In  general,  the  characteristica  of  the  nonlinear  wave  equation  are 
curved  lines  in  the  X-t  plane.  However,  if  a  constant  initial  state  is 
prescribed  for  the  rod.  the  characteristics  of  positive  slope  are  a  family 
of  straight  lines  and  the  con  .'sponding  wave  is  a  simple  (acceleration)  wave. 
If  the  excitation  at  the  end  of  the  rod  is  such  that  successive  wavelets 
are  generated  with  decreasing  shift  rates,  these  straight  characteristics 
diverge  in  the  X-t  plane.  Hut,  if  the  shift  rate  at  the  end  of  the  rod 
increases  (e.g.,  due  to  compression  nr,  sometimes,  sufficient  tension  with 
an  "S-shaped"o  -  X  curve),  the  characteristics  of  positive  slope  will  r.o 
longer  diverge.  Instead,  they  converge  and  form  an  "envelope"  as  shown 
in  Fig.  3.  It  is  on  this  envelope  that  the  values  of  velocity  and  stress. 


carried  by  the  characteristics,  conflict  so  that  the  curve  Cg  is  an 
approximation  of  a  shock  wave  propagating  with  variable  speed  and  carrying 
a  variable  stress.  The  earliest  time  that  such  n  envelope  appears,  i.e., 
wh.-n  the  first  two  characteristics  converge,  a  cusp  is  formed  at  some  point 
XCR.  At  this  point  (XCR,  tCR)  a  unique  solution  of  the  wave  motion, 
characterized  as  a  simple  wave,  is  mathematically  impossible.  It  is  the 
jump  conditions  which  enable  us  to  continue  past  (Xc„,tcR)  with  a  unique 
solution  for  the  shock. 

The  second  category;  the  evolution  of  discontinuities  from  smooth  initial 
data,  is  the  topic  which  seems  to  be  of  current  interest  (see  e.g.,  [17], 
[24],  [25],  [26]).  When  applicable,  the  method  presented  by  Ames  [17,24] 
is  the  simplest  and  the  most  accurate.  This  method  results  from  observing 
that  classes  of  quasilinear  equations  can  be  obtained  by  differentiation  of 
fi-st-order  nonlinear  equations.  The  first  order  equations  are  then  used 
to  calculate  the  time  tCR. 

For  the  one-dimensional  rod  considered  herein,  we  find  from  [17]  that 
the  critical  time  is  given  by 


where  primes  mean  differentiation  of  the  quantity  wtch  respect  to  the  argu¬ 
ment  and  the  forms  of  h  and  depend  upon  which  form  of  the  wave  equation 
wc  are  considering,  (3.1)  or  (3.3).  F>>r  the  displacement  form  of  the  equa¬ 
tion  of  motion,  (3.1),  we  have 

ux  =  h[X  +  c(ux)t]  (3.12) 


0  =  c(ux)  (3. 13) 

with  c(ux)  being  the  material  shift  rate  defined  in  (3.2).  If  we  consider 
(3.3)  on  the  other  hand,  we  have 


/ 
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X  =  h[X  +  0(X)t]  (3.14) 

0  =  c(\)  (3.15) 

where  the  form  of  c(X)  is  given  in  (3.5). 

4.  FINITE  ELEMENT  APPROXIMATIONS 

We  are  now  ready  to  develop  discrete  models  of  the  equations  governing 
nonlinear  waves  by  using  the  finite-element  concept.  Toward  this  end, 
we  begin,  as  usual,  by  partitioning  the  rod  into  a  finite  number  K  of  seg¬ 
ments  connected  together  at  nodal  points  at  their  ends.  The  E+l  nodes  are 
labeled  0  =  X°  <  X1  <  . . .  <  X1 =  L„ ,  and  the  mesh  length  h  »  xa+1  -  Xtt 
(or  =  0,1,...,E)  is  assumed  to  be  uniform.  In  the  finite-element  method, 
we  seek  an  approximate  solution  to  either  (3.1)  or  (3.3)  in  a  finite  dimen- 
alonal  subspace  of  H1(0,L),  the  elements  of  which  are  locally  of  the  form 

f(x,o  -  fa(t>yx)  (4.D 

The  repeated  nodal  index  a  ranges  from  l  to  2  for  the  one-dimensional  rod 
element.  Here,  fQ(t)  is  the  value  of  f(X,t)  at  node  or  of  the  element  at 
time  t,  and  ^(X)  are  the  usual  local  (elemental)  interpolation  functions. 
In  general,  these  local  basis  functions  contain  complete  polynomials  of 
degree  p,  where  p+1  is  the  order  of  the  highest  material  derivative  that 
appears  in  the  energy  equation  for  the  element.  In  this  case,  for  either 
(3.1)  or  (3.3),  we  have  p  =  l  (cf.  [12]). 

To  obtain  approximate  solutions  to  (3.1),  we  take  f(X,t)  in  (4.1)  to 
be  the  local  displacement  field  u(X,t): 

u(X,t)  =  ua(t)1fa(X)  (4.2) 

Approximate  solutions  to  (3.3)  are  similarly  obtained  when  we  take  f(X,t) 
to  be  the  local  extension  ratio  (stretch)  X(X,t): 


/ 
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MX,  t)  =  Xa(t)Vx) 


(4.3) 


4.1  Local  Form  of  the  Equations  of  Motion.  We  consider  first  the  equation 


of  motion  (3.3).  We  seek  approximations  of  weak  solutions  of  (3.3)  by 
requiring  that 


n  n 

/  0dX  =J  0(c2Xx), 


(4.4) 


for  any  arbitrary  function  0(X)  which  has  a  continuous  first  derivative  and 
which  vanishes  outside  the  element.  Integrating  the  right  hand  side  of 
(4.4)  by  parts  gives 


n  n 

/  X«.0<K+/ 


ax  0*dx  =  Lp0  ax  0\ 


(4.5) 


where,  from  (3.2),  da/dX  =  (da/d\)\x  =  c2\x.  If  0*  is  constant,  (4.5) 


becomes 


u 

J  Xtt0dX  +  — 


0*[°3  =  [C?XX0] 


(4.6) 


Taking  a  piecewise  linear  approximation  for  X  in  (4.3),  we  have 


X(X,t)  =  yx)Xff(t)  =  (aff+  baX)Xa 


(4.7) 


where  for  the  ich  element 


laa)  =  [0,1]  ;  {bj  =  -i  [-1,1] 


(4.8) 


Hence,  with  (4.8)  in  (4.7),  we  have 


=  "  Vcr  h(X2  ‘ 

-  b«  ’ -E  t-1'11 


(4.9) 


i 
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Taking  0(X)  =  (X)  in  (A. 6)  and  incorporating  (A. 7),  we  obtain  the 

equation  of  motion  for  a  typical  rod  element 


■VS i + 


<‘•10) 


where  mg^  is  the  consistent  mass  matrix  defined  by 


=  mc*S  *  p»Ao 1  *“*-dx 


u 

*/  W 


(A. 11) 


is  the  generalized  force  at  node  a 


Per  =  A*£*ac*30 

'■*  i'. 


(A. 12) 


Equation  (A. 11)  can  be  integrated  to  get  the  consistent  mass  matrix 
in  the  form 


mPcr  -  6  n(1  +  6p«) 


(A. 13) 


where  m  =  P0A0h  is  the  rasa  of  a  typical  element  and  6„H  i«  the  Kronecker 
delta.  However,  if  the  mass  is  considered  to  be  lumped  at  the  nodes,  ra^ 
will  be  of  the  diagonal  form,  We  prefer  the  lumped  mass  model  in 

this  study,  not  only  for  the  increased  computational  speed,  but  also 
because  it  tends  to  maintain  the  finite  character  of  the  speed  at  which 
waves  are  propagated  along  the  rod  while  simultaneously  reducing  "ringing** 
in  front  of  the  wave  front. 

We  now  turn  to  the  displacement  form  of  the  equation  of  motion  (3.1). 
Analogous  to  the  procedure  used  to  model  (3.3),  we  approximate  the  local 
displacement  field  by  (A. 2) 


u(X,t)  »  *a(X)ua(t)  -  (aa+  bQ)C)ua  (A.1A) 

with  a  and  b  defined  in  (A. 8),  and  here  u  ■  u(X  .t).  Accordingly,  we 
Of  Of  Of  Of 

obtain  as  the  equation  of  motion  for  a  typical  rod  element, 


V#  +  b„A»h3  -  pa 


(4.15) 


In  this  case,  since  ux  =  b^u^,  b°ch  A  ®nd  0  =  a(^)  are  constant  for  each 
element.  The  mass  matrix  is  as  previously  defined,  but  the  generalized 
nodal  force  p  is  now 

Pa  =  A0a*a  (A.  16) 

4.2  Global  Form  of  the  Equations  of  Motion.  To  facilitate  book-keeping, 
we  use  superscripts  as  the  element  index  and  subscripts  as  the  nodal  index. 
If  we  define  P„  as  the  net  generalized  force  applied  at  node  N,  so  that 
PM  =  p!|-1  +  Pi»  then  the  global  equation  of  motion  for  node  N  is  obtained 
from  (4.10)  as 

A0 

PN  =  mXN  -  -jj*  +  Oh+j)  (4,17) 

where,  from  (4.8),  bJJ-1  =  +  h-1  and  bJJ  =  —  h_1x  and  we  have  used  the  lumped 
mass  model.  At  the  ends  of  the  rod,  of  course,  the  form  of  (4.17)  changes 
according  to  the  type  of  boundary  conditions.  Equations  (4.17)  represent 
the  global  system  of  nodal  equations  of  motion  for  the  finite-element 
model  of  the  rod.  Since  the  solution  to  these  equations  will  be  in  terras 
of  the  displacement  gradients  (X  -  1  =  du/dx),  the  nodal  displacements,  if 
desired,  are  obtained  by  spatial  integration. 

The  displacements  can  be  obtained  directly  if  the  global  equations  of 
motion  are  formed  as  above,  but  using  (4.15)  rather  than  (4.10): 

Pn  =  mIl’N  +  MOi-i  .  0N)  (4.18) 

5.  FINITE  ELEMENT/ DIFFERENCE  EQUATIONS 

The  remaining  step  in  discretizing  the  nonlinear  equations  is  approxi¬ 
mating  their  temporal  behavior.  Solving  the  nonlinear  equations  of  motion 
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by  stepwise  numerical  methods  can  be  severely  complicated  by  the  presence 
of  shocks.  Therefore,  we  choose  an  explicit  finite  difference  scheme 
which  automatically  treats  the  shocks,  whenever  and  wherever  they  may 
occur,  without  the  necessity  of  the  tedious  application  of  the  jump  con¬ 
ditions  at  each  time  step  of  the  solution  process.  This  method  is  the 
well-known  Lax-Wendroff  difference  scheme  [27,28],  which  is  generally 
classified  as  an  artificial  dissipative  method.  The  success  of  this 
particular  method  comes  from  applying  finite  difference  approximations 
to  the  governing  equations  expressed  as  ’’conservation"  laws.  (For  details 
of  the  theory  of  conservation  laws  and  weak  solutions,  see,  e.g.,  [23], 

[27],  or  [29].)  The  novel  aspect  of  the  following  temporal  discretization 
is  that,  by  rewriting  (3.1)  as  (3.3),  we  obtain  the  governing  equation  in 
the  form  of  a  conservation  law  which,  with  only  a  linear  finite-element 
approximation,  enables  us  to  develop  a  Lax-Wendroff  type  integration  scheme. 


5.1  Lax-Wendroff-Finite-Element  Scheme.  The  temporal  discretization  is 

accomplished  in  the  spirit  of  the  Lax-Wendroff  equations  (cf.  [28],  p.  302): 
•  •• 

we  first  denote  qH  =  XM  and  expand  qN(t  +  &t)  into  a  Taylor  series  up  to 
second  order  terms 


qN(t+At)  -  qN  (t)  +  AtqN  +  |  (At)aqN  +  0(At3)  (5.1) 

The  t  derivatives  in  (5.1)  are  now  replaced  by  X  derivatives  (except  for 
the  nodal  force  Pt )  by  means  of  (A. 17)  where 


%  “Sh  (°N-1  '  2°N  +  °N  +  1>  +  *  PN 


(5.2) 


"  it  ^  -  a  (V,  - 2°* +  w +  £  K 


(5.3) 


\ 


Since  o  =  o(X),  a  =  P0csq,  so  that  (5.3)  becomes 
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9n  ~  (CN-1^N-1  ”  ^CN^N  +  CN+1^N+1^  +  m  ^N 

Note  that  since  the  quantities  PN  are  prescribed,  so  also  are  the  PN  (e.g., 
if  PN  =  sin  t,  then  PN  =  cos  t).  Substituting  (5.2)  and  (5.4)  into  (5.1), 
we  obtain  the  finite  element/dif ference  equation  for  the  interior  nodes  of 
the  discrete  model: 


c1  -  <x>2  c|  +  t&>a  -  «=;>ax  +  i  3 


AtA 


+  ~3r  <c«*-i  -  2ai +  °;-i> + 15  t2itp; +  <w>ap:J  <5-5> 


where  t  =  nAt,  and  q(t)  =  q(nAt)  =  qn ,  etc.  Similarly,  the  finite  element/ 
difference  equations  for  the  end  nodes  are 


q  r1 


and 


„B+1 

9e+1 


At-  Vi  utn 

(^)aC(ci)%n  +  [fa*  -  (cj)*}q »]  +  2  —  (ol  -  o-) 

+  ^  [2AtP“  +  (At)8p;]  (5.6) 

<£)*l(c't)a q?  +  ((fe)8  -  <C(\x>»qB\x3  +  2  (Of  -  0»+1) 

+  ^2AtPl»+l+  (At)8pf+l]  (5.7) 


To  compute  the  extension  ratios,  we  simply  repeat  the  foregoing  pro¬ 
cedure:  we  first  expand  \£+1: 


*  X*  +  AtqjJ  +  (|^)'q“  +  0(At3) 


and  so,  using  (5.2),  we  get: 


K"  -  N!  +  ««;  +  2p~  <r>aK->  -  2°l +  +  ps 


-  K  +  wq;  +  f  ^  p; 


Kxl  -  +  wq.*..  +  t-  -  °:»>  +  n 


**+i  T  o.  vh 


E  +1  ■ 


m 


E+l 


(5.8) 


(5.9) 
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where  a|J  and  (c” )a  are  determined  from  (3.4)  and  (3.2$>  respectively. 

Equations  (5.5)- (5 .7)  and  (5.9)  are  the  2(E+1)  finite  element/ 
difference  equations  used  herein  for  the  study  of  shock  waves. 

5.2  Velocity-Centered  Central  Difference  Scheme.  To  illustrate  the 
effectiveness  of  the  Lax-Wendroff  method  we  will  compare  the  results  with 
the  finite  element/difference  equations  which  are  obtained  using  the  dis¬ 
placement  equations  (4.18).  The  temporal  discretization  is  accomplished 
for  these  equations  by  using  velocity  formulated  central  difference  [30], 
in  which  the  general  nodal  equation,  u‘N(t)  =  F(t),  is  approximated  by 
introducing  the  nodal  velocity  vN  *  iiN  and  thereby  generating  two  equivalent 
first  order  equations,  which  are  then  differenced  to  obtain 

v"+^  *  +  AtFJj 

"  H  (5.10) 

u”+1  «  uj  +  AtvJ+V 

Direct  substitution  shows  that  (5.10)  is  equivalent  to  using  the  usual 
central  difference  approximation  for  u"N.  However  (5.10)  generally  admits 
to  less  round-off  error  (cf.  [3l]). 


5.3  Stability  and  Convergence  of  the  Central  Difference  Schemes  for  Nonlinear 
Wave  Equations.  We  remark  that  in  a  recent  paper  [3l]  we  presented  a 
detailed  analysis  of  the  numerical  stability  criteria  and  rates  of  con¬ 
vergence  for  finite-element/finite-difference  approximations  of  Lhe  non¬ 
linear  wave  equation  (3.1).  For  completeness,  the  principal  results  of 
that  study  are  summarized  as  follows: 

•  The  stability  of  the  scheme  in  energy  is  assured  if  (h/At)  >  v  C 

Of  max 

v  C  (u.)/2,  where  h  is  the  minimum  mesh  length  for  the  finite-element 

model,  v  (or  =  1,2)  are  constants,  corresponding  to  a  consistent  mass 

formulation  and  v.  to  a  lumped  mass  formulation,  and  C  (u_)  is  the  maximum 
*  max  * 
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speed  of  propagation  of  acceleration  waves  relative  to  the  material. 
Obviously,  this  result  reduces  to  similar  criterion  obtained  for  linear 
difference  approximations  when  CMJ[  =  constant. 

•  To  maintain  stability  for  a  given  fi nite-element/dif ference  scheme 
witii  fixed  h,  it  is  necessary  to  use  a  smaller  time  step  for  the  consistent 
mass  formulation  than  for  the  lumped  mass  formulation  since  v  >  v2 ! 

•  Under  the  stated  assumptions,  the  square  of  the  Lg-norm,  ||e^||2,  of 
the  gradient  of  the  error  at  each  time  step  i  is  0(h2  +  (At)s).  (Similar 
accuracies  are  obtained  after  R  time  steps  in  the  linear  case).  Uniform 
convergence  of  the  error  e  is  also  obtained. 

•  The  same  rates-of-convergence  for  the  consistent  mass  formulation 
are  obtained  for  the  lumped  mass  formulation. 

All  of  these  results  hold  only  at  points  at  which  the  solution  Is 
smooth  -  i.e.,  at  points  not  at  the  wave  front. 

6.  NUMERICAL  RESULTS 

In  thii  section,  we  cite  numerical  results  obtained  from  application 
of  the  preceding  theory  to  representative  problems.  For  our  numerical 
examples,  we  consider  a  thin  rod  of  Mooney  material  (Cv  ■  24.0  psi,  C2  * 

1.5  psi)  with  the  following  undeformed  characteristics:  length  »  3.0  inches, 
cross-sec  >Monal  area  =  0.0314  in2,  mass  density  *  10-4  lb.sec2/in4.  For  the 
finite  element  model,  we  take  60  evenly  spaced  elements,  so  that  Ii  *  0.05 

in.  and  N  =61. 

6 

6.1  Tensile  Loading  (square  wave).  We  consider  a  force  of  constant 
magnitude  applied  a."  the  free  end  of  the  rod  as  a  step  function  at  t  *  0, 


/ 


then  similarly  removed  at  a  later  time  t  =  t*,  i.e.,  a  square  wave.  For 
this  example,  t*  =  0.002  seconds.  Figure  4  shows  the  stress  wave  response 
to  a  two  pound  step  loading  (this  corresponds  to  86%  strain  statically) 
for  both  mass  distributions  and  both  time  integration  schemes  previously 
discussed.  Several  important  items  mentioned  earlier  can  be  observed  in 
Fig.  4. 

(1)  The  lumped  mass  model  tends  to  preserve  the  finiteness  of  the 
speed  of  propagation.  (Note  the  "ringing"  in  front  of  the  consistent  mass 
stress  wave.) 

(2)  The  acceleration  wave  front  does  tend  to  flatten  with  time. 

(3)  The  wave  is  propagated  into  the  undisturbed  portion  of  the  rod 
as  a  simple  wave.  Recall  that  for  a  simple  wave,  Y(t)  la  constant,  so 
that  by  multiplying  Y(0.001)  by  2.0  -  the  ratio  of  the  elapsed  time  incre¬ 
ments  -  we  obtain  Y(0.002),  except,  of  course,  for  that  portion  of  the 
wave  affected  by  the  fixed  boundary. 

(4)  The  Lax-Wendroff  scheme  Is  clearly  superior  to  the  central 
difference  scheme,  particularly  In  the  presence  of  shocks.  Not  only  does 
it  produce  no  ringing  in  front  of  the  wave,  but  the  unloading  shock  wave 

is  depicted  without  the  large  oscillations  behind  the  shock.  (These  oscilla¬ 
tions  have  been  interpreted  by  some  as  numerical  instability  of  the  Integration 
scheme.  This  is  not  so;  the  amplitude  of  these  oscillations  does  not 
change  with  time.  As  pointed  out  in  [28],  these  are  lumped  mass  oscilla¬ 
tions  resulting  from  discretization  error  and  they  represent  the  Internal 
energy  which  must  appear  in  the  "shocked"  region  according  to  the  jump 
conditions.  It  is  conjectured  that  the  Lax-Wendroff  scheme  converts  this 
oscillatory  energy  into  true  Internal  energy.) 

Figure  5  shows  in  some  detail  the  response  of  the  rod  to  this  2- 1 1> 
step  load.  The  results  confirm  the  fact  that  weak  shocks  propagate  in  a 


locity  central  of  fere 
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manner  similar  to  simple  waves  -  the  stress  increases  at  the  wall  almost 
by  a  factor  of  2.0  and  the  stress  wave  is  reflected  from  the  wall  without 
appreciable  change  in  shape  or  magnitude. 

6.2  Sinusoidal  Forcing  Function.  This  example  dramatically  illustrates 
that  the  central  difference  scheme,  without  modification,  cannot  handle 
shocks.  Here  a  concentrated,  time- dependent  load  which  varies  sinusoidally 
is  applied  at  the  free  end;  a  complete  loading  cycle  occurs  in  0.002  seconds. 
It  is  clear  from  the  computed  response  shown  in  Fig.  6  that  shocks  develop 
quickly  for  this  kind  of  loading.  Unlike  the  response  for  the  tensile 
step  load  where  the  unloading  wave  is  produced  by  simply  removing  the 
load,  the  sinusoidal  load  actually  "pushes"  the  end  of  the  rod.  The 
instant  the  load  starts  to  decrease  is  the  moment  when  the  first  wavelet 
is  generated  which  propagates  faster  than  the  preceding  one.  Thus,  at 
some  time  subsequent  to  when  the  compression  cycle  starts,  a  compression 
shock  forms  in  the  rod. 

A  comparison  between  the  two  integration  schemes  is  also  shown  in  Fig. 

6  for  the  sinusoidal  loading.  In  this  case,  it  is  clear  that  the  "shocked 
internal  energy"  behind  the  compression  shock  renders  the  central  difference 
scheme  unacceptable.  It  is  Interesting  to  note,  however,  that  the  tension 
cycle  evidently  "absorbs"  the  large  oscillations  preceding  it  and  again 
produces  a  smooth  wave  front.  The  detailed  response  to  this  loading  is 
shown  in  Fig.  7.  From  the  response  shown,  we  notice  several  interesting 
features  of  nonlinear  wave  motion; 

•  The  compressive  shock  wave  is  reflected  from  the  wall  as  a  compressive 
shock  wave  by  almost  doubling  the  compressive  stress;  but  the  tension 
part  of  the  stress  wave  'a  reflected  with  only  a  small  increase  in  stress. 


Stress 


history  of  stress  wave  response  to  sinusoidal  end  load. 


Fig.  7c.  Time  history  of  stress  wave  response  to  sinusoidal  end  load 


*  At  t  *  A. 7  milllsec,  two  compression  shocks  collide.  The  numerical 
results  shown  here  indicate  that  when  two  shocks  collide  in  a  solid 
material,  they  penetrate  one  another  with  little  or  no  deterioration. 

This  is  apparently  contrary  to  the  collision  of  shocks  in  gases  [21]. 

*  By  comparing  the  response  at  t  =  3  msec  with  that  at  t  =  5  tr.3ec, 
we  note  that  the  response  tends  to  repeat  itself  (with  some  variation 
due  to  the  reflection)  with  essentially  the  same  period  as  that  of  the 
forcing  function. 

*  As  in  the  development  of  shocks  from  Lipschitz  continuous  data,  the 
shock  forms  subsequent  to  iAitiation  of  the  compressive  cycle.  Thus  we  are 
led  to  examine  the  positive  slope  characteristics  in  the  X-t  plane  to  see 
if  they  predict  t„R  for  this  type  of  loading.  Figures  8  and  9  show  thit 
if  we  assume  straight  compression  characteristics  of  positive  slope,  the 
cusp  of  the  corresponding  envelope  in  the  X-t  plane  does,  in  fact,  give 

a  good  estimate  of  the  tc„  observed  in  the  stress -time  plots. 
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